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ABSTRACT 

Context. Observations of edge-on galaxies allow us to investigate the vertical extent and properties of dust, gas and stellar distributions. 
NGC 891 has been studied for decades and represents one of the best studied cases of an edge-on galaxy. 

Aims. We use deep PACS data together with IRAC, MIPS and SPIRE data to study the vertical extent of dust emission around NGC 
891. We also test the presence of a more extended, thick dust component. 

Methods. By performing a convolution of an intrinsic vertical profile emission with each instrument PSF and comparing it with 
observations we derived the scaleheight of a thin and thick dust disc component. 

Results. For all wavelengths considered the emission is best fit with the sum of a thin and a thick dust component. The scaleheight 
of both dust components shows a gradient passing from 70 pm to 250 pm. This could be due to a drop in dust heating (and thus dust 
temperature) with the distance from the plane, or to a sizable contribution (~ 15 - 80%) of an unresolved thin disc of hotter dust to the 
observed surface brightness at shorter wavelengths. The scaleheight of the thick dust component, using observations from 70 pm to 
250 pm has been estimated to be (1.44 ±0.12) kpc, consistent with previous estimates (extinction and scattering in optical bands and 
MIR emission). The amount of dust mass at distances larger than ~ 2 kpc from the midplane represents 2 - 3.3% of the total galactic 
dust mass and the relative abundance of small grains with respect to large grains is almost halved comparing to that in the midplane. 
Conclusions. The paucity of small grains high above the midplane might indicate that dust is hit by interstellar shocks or galactic 
fountains and entrained together with gas. The halo dust component is likely to be embedded in an atomic / molecular gas and heated 
by a thick stellar disc. 

Key words. Galaxies: structure - Galaxies: individual: NGC 891 - Galaxies: spiral - Dust, extinction - Infrared: galaxies - 
Submillimeter: galaxies 


1. Introduction 

Nearby edge-on spiral galaxies offer a special opportunity to 
learn about the vertical structure and properties of dust, gas and 
stellar distributions that extend out from the galaxy midplane. 
Several galaxies of this class have been observed in optical bands 
and show many common properties. The high inclination of 
these galaxies makes them particularly well-suited to conduct 
multi-wavelength studies of the vertical structure of their disc. 
Fitting optical images with radiative transfer models allows in¬ 
formation on the dust and stellar distributions to be obtained. 
Stars and dust are found to follow an exponential distribution 
both radially and vertically: the dust radial scalelength is ~ 1.4 
times larger than that of the stellar distribution while its vertical 
scaleheight is about half of that of stars (jKylafis & Bahcall 1987 


Xilouris et al.|1997||1998||1999||Alton et al.|20041|Bianchi 2007 


Baes et al.|2010[ De Looze et al.|2012 De Geyter et al.|2dT4 1. 

Extending studies further out to galactic haloes, there is in¬ 
creasing evidence for the presence of dusty clouds up to a few 
kpc from the galactic disc that seem to be related to the disc- 
halo cycle (e.g. |Howk|2012| l. Direct optical imaging of edge-on 
galaxies often shows extraplanar filaments of dense clouds back¬ 
lit from stars ( Howk & Savage||1997| 1999| 2000 Rossa et al. 


|2004[ [Thompson et al.||2004[ |Howk||2005| r An estimate of the 

mass of dust contained in these clouds indicate that a relevant 


* Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and with im¬ 
portant participation from NASA. 


fraction of dust can be lifted up to these vertical distances with¬ 
out being destroyed ( |Howk|2005| ). 


Moving even further from the galaxy main body, Zaritsky 


(1994 1 observed dust extinction up to 200 kpc out in the galactic 
haloes of two nearby galaxies, NGC 2835 and NGC 3521. This 
result was supported by Menard et al. P010| l, who measured the 
extinction of background quasars correlating to local galaxies 
and found evidence for dust in the intergalactic medium (IGM) 
up to 1 Mpc from galactic centres. 


GALEX and Swift observations in UV bands by Hodges 


Kluck & Bregman (2014 1 reveal diffuse UV light around late- 
type galaxies out to 5-20 kpc from the galactic disc. The emis¬ 
sion is rather blue and not consistent with light from extraplanar 
stars; the favoured hypothesis is that this emission escapes from 
the galactic disc and scatters off dust in the halo. This finding is 
consistent with the extinction measurements from IMenard et all 
( 2010| l, therefore corroborating the evidence for dust in galactic 
haloes. 

With the advent of Spitzer and more recently of Herschel 
( jPilbratt et al.|201()| ) we have the necessary sensitivity to directly 
observe extraplanar dust thermal emission. Aromatic bands have 


been observed as far as 6 kpc above the disc of NGC 5907 (Irwin 
|& Madden]|2006|l and in galactic winds around several nearhy 


sources ( McCormick et al.|2013| l. Continuum dust emission up 
to a similar height was found in NGC 891 by [Burgdorf et al. 
( |2007l l. 

In this paper we focus on the vertical dust distribution in the 
well-known edge-on galaxy NGC 891. This galaxy is at a dis- 
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tance of 9.6 Mpc (e.g. |Strickland et al.||2004 1 from us and it 
has an inclination of ^ 89.7° ( |Xilouris k al.| 1999| l. The dust 
distribution follows an exponential profile with a scalelength of 
/id ~ 8 kpc and a scaleheight of Zd ~ 0.3 kpc, thinner than the 
stellar distribution but more radially extended (Xilouris et al. 
1999|l . Radiative transfer models by Bianchi ( 2008|l anH Popescu 
et al.| pOOOal |2011[ ), starting from models by [Xilouris et al. 
( 1999| l, are able to reproduce the dust SED; thus the radiation 
field across the galaxy is reliably modelled. The galactic stellar 
distribution comprises the two classical components, an old stel¬ 
lar bulge and an old stellar disc, with the addition of a clumpy 
component confined in a very thin region at the centre of the 
main disc. Bianchi & Xilouris (20111, using the results from 
the 3D radiative transfer model from |Bianchi| ( |2008[ ) are able 
to reproduce most of the characteristics of far-infrared (FIR) 
and submm images obtained with the SPIRE instrument onboard 
Herschel. 

Optical broad-band observations obtained by Howk & 


Savage ( I997|l with the WIYN 3.5-m telescope revealed the 


presence of a complex dusty structure above the galactic plane 
up to a distance z < 2 kpc. Dust emission was later directly de¬ 
tected by Burgdorf et al.| ( 2007| l with the Infrared Spectrograph 
aboard Spitzer at 16 and 22 pm. In particular, the surface bright¬ 
ness at 22 pm drops following an exponential decay with a scale- 
height of 1.3 + 0.3 kpc, consistently with the optical measure¬ 
ments by Howk & Savage (1997 1 . Furthermore, Kamphuis et al. 
( 2007| l using MIPS 24 pm observations detected dust emis¬ 
sion up to ~ 2 kpc from the galactic plane and [Whaley et al. 
( 2009| l, using multi-wavelength observations, found extraplanar 
dust emission and aromatic features up to a vertical distance of 
z > 2.5 kpc. However, Kamphuis et al.| (2007 1 and Whaley et al. 
( |2009| l did not consider in detail the effects of the instruments 
PSF and therefore their measurements are to be considered qual¬ 
itative. 

More recently, [Seon & Witt ( 2012} reported on the discovery 
of a UV halo due to the diffuse dust present above the galactic 
plane of NGC 891 and jSeon et al. ( |2014p , using a 3D radiative 
transfer model, found that the UV halo is well reproduced by a 
model with two exponential dust discs, one with a scaleheight of 
0.2 - 0.25 kpc and a wider component with a scaleheight of 1.2 
- 2.0 kpc. 

In this work, we use Spitzer and Herschel data (from 8 to 
250 pm) to study the extent of dust emission above the galactic 
plane, taking a particular care in the analysis of the effects of 
the instruments Point Spread Function (PSF). We then compare 
the derived dust distribution to the extinction and emission mea¬ 
surements from the literature. This paper is organised as follows: 
Section [^presents the observations and the data reduction used 
to retrieve the data; Sectionj^shows the dust emission profile as 
measured by Spitzer and Herschel instruments. In Section]^ we 
highlight the variability of the vertical profile along the major 
axis of the galaxy; Section shows measurements of the dust 
Spectral Energy Distribution, temperature and optical depth at 
different vertical distances. In Section |6] we discuss our results 
and in Section|7]we draw our conclusions. 

2. Observations and data reduction 


Herschel PACS ( jPoglitsch et al. |2010| l photometric obser¬ 
vations were taken at 70, 100 and 160 pm as part of the 
Herschel EDGE-on galaxy Survey (HEDGES; PI: E. Murphy). 
70 and 100 pm data consist of two cross-scans at a 20" 
s-' scan speed (Ohs IDs: 1342261790, 1342261792 for 70 
pm and Ohs IDs: 1342261791, 1342261793 for 100 pm). 


while 160 pm observations consist of four cross-scans at 
20" s-i scan speed (Ohs IDs: 1342261790, 1342261791, 
1342261792, 1342261793). The//ersc/ze/Interactive Processing 
Environment (HIPE, v. 12.1.0; jOttj |2010| was first used to 
bring the raw Level-0 data to Level-1 using the PACS cali¬ 
bration tree PACS_CAL_65_0. Maps were then produced us¬ 
ing Scanamorphos (v.24.0, Roussel 2013). Images were pro¬ 
duced with pixel sizes of T.'5, 2" and 3", respectively (Fig. [^l. 
We compared these observations with the less deep observa¬ 
tions obtained as part of the Guaranteed Time Key Project 
Very Nearby Galaxy Survey (VNGS; KPGT_cwilso01_l; PI: C. 
D. Wilson) finding agreement between the two datasets. Also, 
deeper observations are available from the open-time program 
OTl_sveilleuJ2 (PI: S. Veilleux). However, the field of view is 
too narrow and border effects would be introduced in the verti¬ 
cal profile; we therefore decided not to use this dataset. 

Herschel SPIRE ( Griffin et al.||2010) l photometric observa¬ 
tions at 250 pm were obtained as part of the VNGS (Ohs ID: 
1342189430). The galaxy was observed in large map mode, cov¬ 
ering an area of 20'x20' centred on the object with two cross¬ 
scans, using a 30 s ' scan rate. Data were reduced with HIPE 
using the SPIRE calibration tree SPIRE_CAL_12_3 and the stan¬ 
dard pipeline destriper to remove baselines. The map was then 
produced using the naive mapmaking procedure within HIPE. 
The resulting image has a uniform background and pixel size 
of 6" (Fig. [^. In-beam flux densities are converted to sur¬ 
face brightness assuming an effective beam solid angle Q. - 
465 arcsec^ (SPIRE Handbook, 201^. SPIRE observations at 
350 and 500 pm were also obtained as part of the VNGS but not 
used here since the instrument resolution at these wavelengths is 
too low for our purposes. 

Spitzer IRAC individual Basic Calibrated Data (BCD) 
frames were taken from the combined observations (8 
Astronomical Observing Requests) acquired by G. Fazio 
(Spitzer Cycle 1, PID 3), and processed with version 18.25.0 
of the SSC pipeline. The data were acquired in High Dynamic 
Range (HDR) mode, so that each frame consisted of one short 
integration and two long integrations. The 442 long-integration 
time BCDs were combined into a single mosaic with MOPEX 
( jMakovoz & Marleau|2005| l. Bad pixels, i.e. pixels masked in the 
data collection event (DCE) status files and in the static masks 
(pmasks), were ignored. Additional inconsistent pixels were re¬ 
moved by means of the MOPEX outlier rejection algorithms. We 
relied on both the box and the dual outlier techniques, together 
with the multiframe reject algorithm. The frames were corrected 
for geometrical distortion and projected on to a north-east coor¬ 
dinate system with pixel sizes of (X'h, roughly half the size of 
the original pixels. The final mosaics were obtained with stan¬ 
dard linear interpolation. The same was done for the uncertainty 
images, i.e. the maps for the standard deviations of the data 
frames. After experimenting with various HDR correction algo¬ 
rithms, we neglected this correction as it resulted in inferior fi¬ 
nal mosaics. Despite several image reduction attempts, the IRAC 
Channel 4 image (see Fig.[2l still shows a residual column bleed, 
visible toward the southwest part of the image. However, these 
artifacts do not affect our analysis as they are masked out when 
necessary. In order to obtain the “dust map” at 8 pm we subtract 
the stellar component following the prescription by jHelou et al.j 
(2004 1 : we scale the IRAC 3.6 pm map by 0.232 and we subtract 
it from the observed IRAC 8 pm map on a pixel-by-pixel basis. 


* http://herschel.esac.esa.int/Docs/SPIRE/pdf/spire_ 
om.pdf 
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Fig. 1: From top to bottom: IRAC 8 jum (after stellar subtraction), MIPS 24 pm, PACS 70 pm (left column) and PACS 100 and 160 
pm and SPIRE 250 pm (right column) images of NGC 891. Each panel shows an area of ~ 14' x 7' centred on the galaxy’s radio 
continuum coordinates a = 2^22™33'®.0, d = 42°20'57'.'2 (J2000.0, blue cross; Oosteiioo et al. 2007). The main beam size is indicated 
by the filled blue circle (EWHM: 2", 6", 6", 7.'2, 11", 17.'8) and images have pixels sizes of a.'6,l'.'5, 1'.'5, 2", 3", 6". The white 
arrow at the top-right corner shows the Z-axis (i.e. the axis parallel to the short side of the instrument field of view) of the spacecraft 
during observations. The sky RMS noise is cr = 0.03,0.03,1.06,1.23,0.78 and 1.31 MJy sr '. Contours are shown at 3, 10, 100, 
200 and 1000-cr over all images. North and West directions are indicated. 


Spitzer MIPS data at 24 pm were processed by Bendo et al. 
(2012) and the final image has a pixel size of T.'5 (Pig.[2l- 

PACS and SPIRE images were created on a grid rotated of 
an angle equal to the galaxy’s position angle (PA = 22.9, as esti¬ 
mated by PA fitting performed by Bianchi & Xilouris|2011| and 
Hughes et al. 2014|l in order to have the galactic disc on the hor¬ 


izontal direction. IRAC and MIPS data were rotated by the same 
angle in post-processing after data reduction. Eor each image, 
the mean sky flux estimated from different background apertures 
was subtracted. 


consider a general vertical profile of the form: 

7d(^) = 7d,i exp(—^ I -H /d ,2 exp(— —] , (1) 

\ ^d,l/ \ 2d,2/ 

where i -i- /d ,2 is the dust surface brightness at the galactic mid¬ 
plane, and Zd,i and Zd ,2 are the scaleheights of two dust compo¬ 
nents. We then convolved this general vertical profile with the 
PSP of each instrument used in our analysis and compared the 
resulting profiles with observations. 


3. Dust emission profiie 

In order to test the possibility to quantify emission both from 
a thin dusty disc and from a thicker halo dust component we 


3.1. Point and “Line" Spread Functions 

For an extended object, the observed surface brightness distribu¬ 
tion, g(i, ]), is the result of the convolution of an intrinsic surface 
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brightness distribution, /(/, j), with the instrument PSF(!, j), 
g(/,;) = [PSF */](/,;•) 

Nil! Nj/2 

^ ( 2 ) 

s=-Ni/2 t=-Njl2 

where i and j are two perpendicular directions and Ni and Nj are 
the number of pixels along these two directions, respectively. A 
(resolved) edge-on galaxy is a particular case of extended emis¬ 
sion source where the surface brightness has a smooth gradient 
along the galactic disc compared to the perpendicular direction 
(here defined as i and j directions, respectively). Assuming that 
the intrinsic emission is constant along i and infinitely extended 
(or at least much more extended than the considered region), so 
that f(i, j) = f{l, j) = f(j), with 7 any position along 7, we can 
write: 


Nj/2 


Ni/2 


8(1 j)- X Z PSF(7-.,;-0 

t=-Njl2 L5=-MV2 

Njl2 

- 2 /(f)LSF(;-f) 

r=-Nj/2 

where the “Line ” Spread Function (LSF) 


LSF(;) = £ PSF(7,;) 


( 3 ) 


represents a vertical cut of the two-dimensional image resulting 
from the convolution of the instrument PSF with a infinitely thin 
line. Under the above assumptions (see Appendix for a test 
on the assumed geometrical model), a vertical profile extracted 
from the image at a position 7 can be compared equivalently to 


- a profile, g(J, j), extracted from a modelled image obtained 
from the convolution of the two-dimensional intrinsic sur¬ 
face brightness distribution, f(i, j), with the PSF (Eq.|^; 

- a one-dimensional model, f{j), convolved with the LSF 
(Eq.§. 


The first appr oach was used by [Bianchi & Xilouris|(|201 1 , while 
the second by Verstappen et al. ( |20131 l arid Hughes et al. ( |2014| l 
(though without a formal demonstration of its validity). In this 
work we use the latter. 


3.2. Spitzer and Herschel PSFs 

Knowing the PSE of an instrument with a good precision is key 
to perform convolution (or deconvolution) operations on images 
and we therefore carefully consider several details to have the 
best estimate of each PSE. In particular, the vertical profile of 
this source is resolved but its extent is comparable to the EWHM 
of the PSE, therefore making the choice of the right set of PSEs 
a critical issue. 

Theoretical Herschel PACS PSPs are released by the PACS 
Instrument Control Centre (ICC) assuming a spectral slope a - 
-4 in F.i, which corresponds to the ideal PSP that would have 
a point source with a very high temperature that emits with 
a blackbody spectrum, By(T) (we will refer to these PSPs as 
“A(std)”). However, dust emits with a modified blackbody spec¬ 
trum, (v'f ByiT), where p is the spectral index and By{T) is a 
blackbody radiation for an equilibrium temperature T. Dust in 
galaxies is expected to be in thermal equilibrium at an average 


temperature T ~ 20 K and a spectral index p ~ 2, leading to 
a modification of the effective PSP. We therefore modified the 
A(std) PSPs according to the dust spectrum assuming T = 20 K 
and P - 2 (see Appendix |B.1| for a detailed description of the 
method used). We will refer to these PSPs as “A(mod)”. 


We compared both the standard and modified PACS PSPs 
with observations of Vesta and Mars (point-like sources in PACS 
bands) and we noted relevant differences in the emission profile. 
Por this reason, we combined these observations to have a better 
estimate of the observed PSP (“B(std)”, see Appendix B.2 for 
details). As for the theoretical PSPs from the ICC we modified 
the observed PSPs according to the dust spectrum assuming the 
same temperature and p values (“B(mod)”). 


Herschel SPIRE PSPs are released by the ICC and are ob¬ 
tained from scan-map data of Neptune and represent therefore 
empirical beam maps, “A(std)”. Pollowing the same method as 
for the PACS PSPs, we modified them according to the dust 
spectrum assuming the same temperature and p values (see 
Appendix |B.l I for details, “A(mod)”). 

The IRAC Instrument Support Team noted that IRAC PSPs 
are undersampled and then developed Point Response Punctions 
(PRPs) which combine the PSP information and the detector 
sampling. The PRPs were generated from models refined with 
in-flight calibration test data involving a bright calibration star 
observed at several epochs, making them good estimates of the 
observed bearr0 


Th eoretical PSPs for th e Spitzer MIPS camera were gener¬ 


Aniano et al. 


(2011 1 using the software STinyTinj^They 


ated by 

assumed a blackbody source at T = 25 K and used a pixel size of 
0.5" (we will refer to the MIPS 24 pm PSP as “A(std)”). As the 
optics are very smooth, the PSPs generated in this way matches 
the observed PSPs with a sufficient precision (MIPS Instrument 
Handbook). 


After several tests, we define our best estimate of PSPs as 
A(std) for IRAC and MIPS, B(mod) for PACS and A(mod) for 
SPIRE. All these PSPs cover a large dynamic range (> 10®) and 
are defined until large radii from the centre (> 3 arcmin); thus 
they should not miss any faint point-source extended emission. 


None of the PSPs described above have been circularised 
and the relative orientation between PSP and image needs to be 
taken into account. This is achieved by rotating the PSP image 
so that the direction of the spacecraft Z-axis aligns with the di¬ 
rection of the spacecraft Z-axis on NGC 891 images (indicated 
by the white arrow for all images in Pig. [^. Subsequently, to 
obtain the LSP, for each vertical distance, pixels are summed 
along the direction of the galactic major axis. The LSP is always 
broader than the PSP because of the integrated contribution of 
the PSP Airy rings off the main beam (PSP and LSP would have 
the same PWHM, instead, if the PSP were Gaussian and circu¬ 
lar). However, some positions along the major axis of an edge- 
on galaxy might be dominated by the emission of unresolved 
sources (such as bright star-forming regions), rather than from a 
more diffuse disc. The vertical profile above those regions would 
then appear narrower. We will discuss the impact of this issue on 
our analysis in Section A vertical cut of PSPs and PSPs for 
all the considered instruments are shown in Pig. in blue and 
orange dotted lines, respectively. 


^ IRAC PRFs are available at http: //irsa. ipac. caltech. edu/ 
data/SPITZER/docs/irac/calibrationfiles/psfprf 
^ MIPS PSF FITS files are available at http://www.astro. 
princeton.edu/~ ganiano/Kernels.html 
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z (kpc) z (kpc) 

Fig. 2: Normalised vertical profiles for observations (black solid lines). The green dashed lines are the result of the convolution of 
the one-component dust emission prohle with the instrument LSF while the red dashed lines are the result of convolution of the 
two-component dust emission profile. Normalisation values are 82.03, 48.17, 1119.47, 2023.22, 1798.00, 630.58 MJy/sr for IRAC 
8, MIPS 24, PACS 70, 100, 160 and SPIRE 250 jim observed prohles, respectively. Fit results are indicated in the green and red 
boxes, respectively (scaleheights expressed in kpc and the intrinsic surface brightness in MJy/sr). LSF and PSF vertical prohles are 
illustrated in blue and orange dotted lines, respectively. 


3.3. Vertical profile extraction and error estimation 

In order to extract vertical prohles from observations we con¬ 
sider strips of (-7,-1-7) kpc parallel to the major axis of the galaxy 
(corresponding to -2.5 ',h- 2.5') from the centre of the galaxy (a 
= 2‘'22”33*.0, 6 = 42°20'57.'2, J2000.0, blue cross in Fig. Q, 
where the radial prohle along the midplane is relatively constant 
(emission always larger than ~ 1/10 of the peak, see Fig.|^ in 
all bands considered. For each distance from the midplane (the 
vertical step is one pixel) we then calculate the median emission 
over all the pixels in the row within the specihed area. Vertical 


prohles are normalised to unity and shown in Fig. in black 
solid lines. 


For a comparison of the convolved dust emission prohle with 
observations a careful estimation of the uncertainties is needed. 
Following Ciesla et al. ( 2012| l, Roussel ( 2013| l and Cortese et al. 


( |2014| l we identihed two main sources of uncertainties. The hrst 
source is the instrumental noise, o-jnstr, based on the error map. 
It depends on the number of scans crossing a pixel and can be 
estimated using the sum in quadrature of the error map pixels in 
the measurement region. The second source is from background 
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noise, cr^i^y, and can be estimated as: 


cr 


2 

sky 


= A^apCr 


2 

skypix 


• ^apO- 


2 

skymean ’ 


(4) 


where o-s^ypix is the RMS of the fluxes in the chosen aperturfl 
A^ap is the number of pixels in the aperture and CTskymean is the 
standard deviation of the mean value of the sky measured in dif¬ 
ferent apertures far from the galaxy. For each vertical distance, 
the aperture is defined as the row of pixels over which we cal¬ 
culate the median emission (in this case from -7 kpc to +1 kpc 
from the galactic centre). 

We then obtain the total uncertainty, crtot, as: 


It has to be noted that, since we are not directly comparing im¬ 
ages at different wavelengths, calibration errors are not consid¬ 
ered here but will be later in Section |5] 


3.4. Fitting resuits 


We consider two possibilities, a single dust component (Ic), for 
which /(j 2 = 0 in Eq. [^or two separated dust components (2c). 
We then convolve the intrinsic vertical profile to the correspond¬ 
ing instrument LSF (we use our best set of PSFs here) and used a 
-minimisation routine to fit the observed vertical profiles. The 
resulting convolved vertical profiles are shown in Fig.j^in green 
(Ic) and red (2c) dashed lines, normalised to unity at the galac¬ 
tic midplane. Fitting parameters are indicated for the two case^ 
the intrinsic dust surface brightness is expressed in MJy/sr and 
scaleheights in kpc. Both visually and from fitting parameters, 
we note that for all considered instruments the two dust compo¬ 
nent case best fit the observations. In particular, = ;^f^/DOF 
for the two dust component case is < 5 (with the exception of 
IRAC 8 fim data where ~ 30), while for the case of a single 
dust component is always > 5 times larger. Using different sets 
of PSFs, the assumption of two dust components always leads to 
a better fit to observations for all the considered wavelengths. 

Hereafter we therefore consider only the case of two dust 
components. We have convolved the intrinsic vertical profile 
with all the PSFs described earlier and we report the scaleheights 
of the thin and thick components for different PSFs and for all 
the considered instruments in Table [T] The errors indicated in 
this table result from the ;^f^-minimisation technique. Region X 
refers to a particular region in NGC 891, at a projected radial dis¬ 
tance hx ~ 5.6 kpc SSW from the galactic centre, where the in¬ 
fluence from unresolved sources was estimated to be lower (see 
Section |^. Scaleheights for this region are illustrated only for 
our best estimates of the PSFs. 

While the use of standard (std) or modified (mod) PSFs does 
not affect the scaleheights, a substantial change is achieved using 


cTsjjypj,; defined in this way includes not only the RMS of the back¬ 
ground sky but also possible variations of the flux due to the non¬ 
flatness of the radial profile of the galaxy at the different distances from 
the midplane. _ 

Hughes et al. 1 2014| also performed single components fits at 100, 


160 and 250 pm, finding vertical scalelengths larger by up to a few 
times those shown in Fig.|^of this work. However, an inspection of the 
vertical profiles presented in Fig. of their paper, and in particular of 
their LSF, suggests that the vertical scales have been accidentally multi¬ 
plied by constant factors. If the factors were removed, their scalelengths 
would be 0.14, 0.15 and 0.23 kpc at 100, 160 and 250 pm, respectively. 
The values are in much better agreement with our estimates, if we con¬ 
sider also the different treatment of the PSF we do in this work. 


A or B PSFs for PACS (see Appendix |B.3| for more details). 
Moreover, scaleheights for region X are always larger (by more 
than a factor two in some cases) than estimates integrated over 
the global (±7 kpc) strip described in Sect. |3.3| 

Fig.j^shows the two scaleheights (black and red symbols for 
thin and thick components, respectively) for our best set of PSFs. 
Error bars indicate the errors resulting from the;^f^-minimisation 
technique. Results from fitting optical/NIR images with radia¬ 
tive transfer models are shown in green triangles (|Xilouris et al. 


19991. Red and black dots are the dust scaleheights for region X. 


Scaleheights obtained from observations at galactic scale 
show a wavelength gradient for A ^ 70 pm for both the thin 
and thick components. However, estimates from region X show 
a smoother gradient for the scaleheight of the thin component 
and an almost constant scaleheight for the thick component 
with variations of < 30%. The observed gradient in the scale- 
height of the thin component may indicate that (see discussion 
in Section]^: 


1. the dust temperature quickly drops within its scaleheight 
leading to the dimming of the FIR emission or 

2. also region X is affected by the emission of point-like 
sources that emit at relatively warm temperatures and there¬ 
fore tend to reduce the effective scaleheight depending on the 
relative wavelength. 

Furthermore, the scaleheights of the thin component ob¬ 
tained from 8 and 24-pm emission are very similar both from 
measurements at galactic scale and in region X. This indicates 
that these estimates are affected by the presence of point sources 
to the same degree. On the other hand, the thick component at 
24 pm is ~ 1.6 times more extended than at 8 pm. 



Fig. 3: Scaleheights of the thin (black diamonds) and thick (red 
squares) dust components. Green triangles represent the dust 
scaleheights as estimated by |Xilouris et al.| ( fl999| l. Black and red 
dots are the scaleheights as measured in region X (see Section]^. 
Our best set of PSFs is used for convolution. 


4. Vertical profile variation 

To inspect the variability of the vertical profile across the galaxy 
we subdivided the midplane into radial regions. Each region has 
been chosen to have a radial width that is twice as large as the 
FWHM of the instrument PSF. Radial profiles for all considered 
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Table 1: Dust emission scaleheights (for thin and thick discs) as measured using IRAC, MIPS, PACS and SPIRE instruments. The 
PSFs used for the convolution are indicated by A/B (std, obs), see text for details. 



IRAC 8 pm 

MIPS 24 pm 

PACS 70 pm 


Zd.l 

Zd,2 

Zd.l 

Zd,2 

^d,l 

Zd,2 


(kpc) 

(kpc) 

(kpc) 

(kpc) 

(kpc) 

(kpc) 

A(std) 

0.114 ±0.01 

0.575 ± 0.001 

0.106 ±0.001 

0.96 ± 0.01 

0.119 ±0.002 

0.76 ± 0.03 

A(mod) 

- 

- 

- 

- 

0.114 ±0.002 

0.76 ± 0.03 

B(std) 

- 

- 

- 

- 

0.074 ± 0.002 

0.56 ± 0.02 

B(mod) 

- 

- 

- 

- 

0.054 ± 0.003 

0.52 ± 0.02 

Region X 

0.145 ±0.01 

0.822 ±0.01 

0.140 ±0.001 

1.36 ±0.01 

0.085 ± 0.007 

1.28 ±0.32 



PACS 100 pm 

PACS 160 pm 

SPIRE 250 pm 


Zd.l 

Zd,2 

Zd.l 

Zd.2 

Zd.l 

Zd.2 


(kpc) 

(kpc) 

(kpc) 

(kpc) 

(kpc) 

(kpc) 

A(std) 

0.138 ±0.002 

0.88 ±0.03 

0.191 ±0.002 

1.17 ±0.02 

0.218 ±0.008 

1.10±0.13 

A(mod) 

0.135 ±0.002 

0.90 ± 0.03 

0.188 ±0.002 

1.21 ±0.02 

0.214 ±0.008 

1.03 ±0.13 

B(std) 

0.101 ±0.002 

0.67 ± 0.02 

0.115 ±0.003 

0.91 ±0.01 

- 

- 

B(mod) 

0.087 ± 0.003 

0.63 ± 0.02 

0.100 ±0.003 

0.89 ±0.01 

- 

- 

Region X 

0.114 ±0.004 

1.02 ±0.10 

0.140 ±0.003 

1.14 ±0.03 

0.239 ± 0.007 

1.40 ±0.12 


wavelengths are shown in Fig. and shaded regions indicate the 
radial portions into which the midplane has been divided. The 
position of the main peaks in emission is highlighted by the blue 
dotted lines. 


We define “region X” the vertical strip cutting the midplane 
at a projected distance hx SSW from the galaxy centre. This 
region has been chosen to have a radial width twice as large 
as the FWHM of the instrument PSF and not to comprise any 
of the main emission peaks for all the wavelengths considered. 
Furthermore, the continuum-subtracted Ha emission-line image 
of this galaxy (Fig. 10 in Howk & Savage|2000 1 shows that most 
of the galactic disc is bright in Ha, indicating the presence of 
several large Hii regions. However, at the position corresponding 
to our region X, the Ha emission is the lowest observed in the 
midplane, therefore suggesting a low contamination by bright 
and hot dust emission from Hii regions and making this region 
unique in this source. In Fig. we indicate region X with green 
dashed lines. 


We extract the vertical profile for each of the considered re¬ 
gions and for all the considered wavelengths. In Fig.|^we show 
the vertical profiles for PACS 100 jim as an example. The in¬ 
tensity of the shades of grey indicates the surface brightness 
of dust emission (black being the brightest region). A vertical 
cut of the PSF and of the LSF are indicated by orange and blue 
lines, respectively. The convolved vertical profile obtained from 
fitting the observed vertical profile of the whole galaxy (-7 kpc 
to 7 kpc from the centre) is shown in red dotted line, while the 
green dashed line indicates the convolved vertical profile ob¬ 
tained from fitting observations in region X. 

Our analysis shows that there is some variation in the verti¬ 
cal profiles. Bright regions have a narrower vertical profile indi¬ 
cating that these regions might be dominated by isolated point 
sources or by a very thin disc not resolved in the vertical profile 
(see Appendix [a] for further details). On the contrary, in region 
X, where the radial profile in the midplane is rather flat, the ex¬ 
tracted vertical profile is close to the widest observed vertical 
profile. Region X is therefore affected to a low degree by the 
emission of unresolved sources and for this reason the vertical 
profile extracted in this region represents a better estimate of the 



-10 -5 0 5 1 

R (kpc) 


Fig. 4: Radial profiles for IRAC 8 jum, MIPS 24 urn, PACS 70, 
100 and 160 yumand SPIRE 250 yum(from top to bottom). The 
position of the main peaks in emission (blue dotted lines) and 
the region X (green dashed lines) are indicated. 


real extent of extra-planar dust with respect to the vertical pro¬ 
file extracted at galactic scale. Scaleheights of the thin and thick 
discs retrieved from fitting observations in region X are indicated 
in Fig. [3] (see the black and red dots, respectively) and listed in 
Table 


5. Dust Spectral Energy Distributions 

By fitting the median vertical profiles of region X for each 
wavelength, we can construct dust Spectral Energy Distributions 
(SED) for each vertical distance from the galaxy midplane. We 
computed the uncertainties on the surface brightness in this way: 
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z (kpc) 

Fig. 5: Vertical profiles at 100 /rm as measured in different radial 
regions. A vertical cut of the PSF (orange) and of the LSF (blue) 
is indicated. The vertical profile for the whole galaxy and for 
region X are shown in red and green, respectively. 


dust SED obtained from fitting the dust model while dotted lines 
represent the fit to a modified blackbody with fixed 13-2. The 
emission intensity decreases with increasing z and the modified 
blackbody peak position does not follow a clear trend. 



wavelength (|xm) 


considering the parameters for region X listed in Table[T] we ran¬ 
domly generated a large number of profiles following a Gaussian 
distribution centred at z^.i and Zi .2 and with the standard devia¬ 
tion given by the error on parameters. For each vertical distance 
we then computed the mean intensity and the standard devia¬ 
tion of the distribution which will represent the error on the 
mean intensity. Furthermore, since we are comparing fluxes in 
different bands, we need to include the calibration errors, cTcai, 
which will be added in quadrature to the errors computed earlier. 
We consider a calibration uncertainty of 3% (IRAC Instrument 
Handbook), 4% ( |Bendo et al.|2012] l, 7% (|Balog et al.|2014|) and 
7% (SPIRE Handbook Version 2.5) for IRAC, MIPS, PACS and 
SPIRE data, respectively. 

5. 1. Dust modelling 

Eor each vertical distance we fit the intrinsic spectrum from 8 
gm to 250 fj.m using the Jones et al. (2013) dust model and the 
DustEM cod^to compute the dust SED. During this procedure 
we use the standard [Mathis et al.| ( |1983] l interstellar radiation 
field and let the following parameters vary: 

- Go, a parameter indicating the intensity of the radiation field, 

- fitmin, the minimum size of small grains in the model, 

- msQInline, the ratio between the mass of small grains and that 
of large grains. 

In this way we can constrain the radiation field intensity and the 
fraction of small grains. 

Eurthermore, in order to have information on the average 
dust temperature and on its opacity, we fit the SED from 70 gm 
to 250 fim to a single-temperature modified blackbody radiation 
(MBB): 

G = Tyg{vlvo'fBy{T), (6) 

where Ty„ is the optical depth at a frequency vq corresponding to 
100 gm, p is the spectral index at vo and By{T) is a blackbody 
radiation for an equilibrium temperature T. 

Fig. 0 shows the spectra at vertical distances z - 
0,0.6,1.5,2.5 kpc (from top to bottom). Solid lines represent the 

® http://WWW.ias.u-psud.fr/DUSTEM/ 


Eig. 6: Spectra for z = 0, 0.6, 1.5, 2.5 kpc (top to bottom). The 
fit using the Jones et al. (2013) dust model is indicated in solid 
lines. Spectra from 70 /urn to 250 /urn are fit with modified black¬ 
body radiation with fixed P - 2 (dotted lines). 


Eurthermore, we find that the value of Amin goes from a^iin = 
(0.40 + 0.05) nm at the galactic disc to Omin - (0.70 + 0.05) nm at 
2.5 kpc above the disc. Consequently, the mass of small grains 
is reduced from mscImi^Q = (6.2 + 0.1) x 10“^ to mscImi^G = 
(3.9 + 0.1) X 10“^. This indicates that small grains are partly de¬ 
stroyed when moving from the galactic disc up to 2.5 kpc above 
it. This could suggest that dust is entrained by interstellar shocks 
or galactic fountains and during coupling to the gas small grains 
are partly eroded because of the high relative velocity with the 
gas. Similar results have been observed in supernova remnants 
in the Large Magellanic Cloud (Borkowski et al.|2006 


et al. 


et al. 


2006|l and obtained from theoretical modelling (Micelotta 


20T0t[Bocchio et al.|2014] i. 


Williams 


As a further analysis, we include dust collisional heating due 
to the presence of the hot halo gas, with an updated version of 
the DustEM code as described by Bocchio et al.|([2013| l. Erom 
deep Chandra data, [Hodges-Kluck & Bregman[(2M3|l found hot 


halo gas at a temperature Tj 
exponential distribution: 


gas 


0.2 keV following a vertical 


n(z) = no exp(-z/zhg) (7) 

where no ~ 6 x 10“^ cm“^ and Zhg = 5 + 2 kpc. Large grains are 
not affected by this extra source of heating and their temperature 
is completely determined by the radiation field intensity. On the 
other hand, small grains reach high temperatures and will emit 
more than if heated only by the radiation field. This leads to 
flmin = (0.90 ± 0.05) nm and mscIniLC - (3-3 ± 0.1) x 10“^ at 
2.5 kpc above the disc, supporting the idea that small grains are 
partly eroded moving from the galactic plane to 2.5 kpc above 
it. 


5.2. Dust temperature estimation 

One of the main parameters that we obtain from fitting data to a 
MBB is the luminosity-weighted mean dust temperature at each 
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Z (kpc) 



Fig. 7: Temperature (left panel) and Tg(/r = hx) (right panel) profiles as obtained from MBB fitting (black lines). Contributions from 
the thin (dotted lines) and thick (dashed lines) discs are shown. Temperature profiles obtained using the Jones et al. (20131 and the 
radiative transfer models RTl, RT2 and RT3 are indicated in blue, red and green lines, respectively. In the inset we show Tg(/t = hx) 
for the thick component and indicate its scaleheight. 


vertical distance which is shown in the left panel of Fig. in 
black solid line. We also show the dust temperature for the thin 
and thick disc separately in dotted and dashed lines, respectively. 
The dust temperature in the thin disc drops off rapidly while the 
dust temperature in the thick disc stays more elevated further 
from the midplane. However, the dust temperature estimation 
depends strongly on the obtained scaleheights for the two dust 
components. In particular, the steep gradient in the scaleheight 
of the thin component might be affected by emission from a thin, 
unresolved, disc and this effect would be reflected on the temper¬ 
ature estimation for z < 1 kpc (see discussion in Sect. 6.1 1 . 

As a comparison, exploiting the results from radiative trans¬ 
fer models and assuming a dust model it is possible to calcu¬ 
late the resulting dust temperature. We assumed the dust optical 
properties and size distributions from Jones et al. ( 2013| l and we 
use different radiation fields, specific for this galaxy, obtained 
from radiative transfer models from the literature: 


Popescu et al.|(|20lT]l;|Popescu & Tu^(|2013 1, RT1 


Bianchi (|2008[ ); Bianchi & Xilouris| ( |2()i l| l, RT2 

Bocchio et al. ( |2012| l, RT3 


Using the DustEM code we then compute the dust SED and fit¬ 
ting it to a MBB for T > 70 jim we estimate the average temper¬ 
ature of the large grains. 

The resulting dust temperature profiles are shown in Eig. 
for RTl (blue line), RT2 (red line) and RT3 (green line). It should 
be noted that, using the Draine & Li ( 2007| ) dust model (which 
was originally used in RTl and RT2) instead of that by |Jones 


|et al.| ( |20T3| ) the temperature is systematically higher but the dif¬ 
ference between the two estimates is never higher than ~ 1 K 
and that temperature profiles follow the same trends. 

The dust temperature profile given by our fit follows the 
same trend as that estimated using radiative transfer models only 
for z > 1.5 kpc where the dust temperature is dominated by the 
thick dust component. On the contrary, the dust temperature of 
the thin dust component drops off more quickly than expected, 
reaching T ~ 16K at ~ 0.5 kpc from the galactic plane. In 
the case region X is not affected by the emission of unresolved 
sources, this discrepancy would require the presence of two dis¬ 
tinct heating sources for the two dust component (see Sect.|6.2|i. 


5.3. Optical depth estimation 

The other main parameter obtained from the fitting procedure 
is the optical depth at 100 pm at projected radial distance hx, 
= hx). In order to compare this parameter with the face-on 
central optical depth in the B-band, from the three radiative 
transfer models cited above (i.e. RTl, RT2 and RT3) w e first 
convert Tinn in tb using the - 7.5 x 10 

et al. (2013J dust model. A similar ratio, tioo/tb 


from Jones 

__ _ 1.0 xlF^ 

IS assumed by the |Draine & Li|p007[ ) dust model. The profile of 
Tg(/i = hx) is shown in the right panel of Eig. in black lines 
and contributions from the thin and thick dust components are 
indicated with dotted and dashed lines, respectively. A zoom of 
the thick component is illustrated in the inset. 

Then, following Kylafis & Bahcalf] ( 1987| l, the conversion 
between the T^(h = hx) and is given by: 

T^(/7 =/lx) = 18.97 (8) 

We see from Eig. |7]that Tg(/! = hx) - 40 at the midplane. Using 
the conversion factor from Eq. |^we then obtain = 2.1. The 
typical values of for radiative transfer models of this galaxy 
are about a factor two higher: = 

and 


4 for 


(|^ 


Bianchi & Xilouris 


Seon et al. 


(20111 and - 3.5 for 


(20141 


Popescu et al. 


We have also performed a fit of the optical depth profile to 
an exponential function of the form: 


\{h = hx,z) = Tg(/! = hx,z = 0) exp(-z/Zd, 2 ) 


(9) 


This fit gives a scaleheight of Zd.i = (1.44 + 0.12) kpc, corre¬ 
sponding to the thick dust component, in good agreement both 
with the estimate by Burgdorf et al. (20071, Zd ,2 - (1.3 +0.3) kpc 
and with that by |Seon et al. j 20T4| l, Zdi ~ [1-2 - 2.0] kpc. Also, 
we calculated that ~ 2 - 3.3% of the mass of dust is present fur¬ 
ther than 2 kpc from the midplane, which is consistent with the 
estimates by Seon et al. ( |2014[ ). 

Einally, we compare the value of Ay as a function of the 
vertical distance with the average value observed by [Menard 
jet al.j ( |2010] l from halo extinction measurements towards dis¬ 
tant quasars. Erom our observations we find Ay = 0.28,0.13 
and 0.064 at z = 1,2 and 3 kpc, respectively, in good agreement 
with values extrapolated from results by Menard et al. (2010): 
Ay = 0.20,0.11 and 0.076 at these three distances. 


9 













































































































































M. Bocchio, S. Bianchi, L. K. Hunt, R. Schneider: Halo dust detection around NGC 891 


6. Discussion 


As seen in Section 3.4 the data for NGC 891 require the pres¬ 
ence of two dust components; this result does not depend on the 
choice of the PSF used in the convolution process. Furthermore, 
in order to limit the contamination from isolated point sources 
or from an unresolved disc, we chose to focus on region X 
for the estimation of the vertical profile. Fig. shows that the 
scaleheight of the thick dust disc in region X is almost con¬ 
stant and undergoes wavelength variations of < 30%, while the 
scaleheight of the thin disc is always smaller than estimates by 
Xilouris et al.| ( T999| and presents a gradient in wavelength and 
thickens by a factor of ~ 4 from 70 to 250 /im. A similar 
trend is seen for all the considered PSFs (see Appendix |B.3| l and 
it is therefore independent on this choice. 


6.1. Thin dust component 


The observed gradient in the scaleheight of the thin component 
has strong consequences on the dust temperature estimation. The 
dust temperature of the thin component drops off very quickly 
within the first 0.5 kpc from the midplane making the optical 
depth estimation more difficult (see for example the fit to the 
SED at z = 0.6 kpc in Fig. |^. Furthermore, the retrieved dust 
temperature for z ~ 0.1 - 0.8 kpc is significantly lower than the 
dust temperature estimated from radiative transfer models. 

Although our estimates of scaleheights are based on obser¬ 
vations of region X, which is the region least affected by the 
presence of point sources or of an unresolved disc, even this re¬ 
gion might suffer contamination from the emission of unresolved 
sources. This would result in a narrower vertical profile and this 
effect is likely to be wavelength-dependent, possibly leading to 
a gradient in scaleheights like that observed. 

To test this hypothesis, we add to the standard two- 
component intrinsic vertical profile an unresolved super-thin 
disc with a scaleheight fixed to be equal to half of the observed 
image pixel size (13-45 pc). We then convolve the intrinsic ver¬ 
tical profile to the instrument LSF and fit the resulting profile to 
the observed vertical profile in region X. 


Our fits show that the relative contribution of the super-thin 
disc with respect to the thin and thick components is wavelength- 
dependent and becomes less important going from 70 pm to 
250 pm; the relative contribution of the super-thin disc to the 
observed midplane emission ranges from ~ 80% at 70 pm to 
~ 15% at 250 pm. This indicates that the spectrum of the un¬ 
resolved emission from the midplane is not grey but it peaks 
at wavelengths < 70 pm, corresponding to a dust temperature 
> 40 K, consistent with the expected temperature range of dust 
heated by nearby hot stars. The super-thin disc could thus repre¬ 
sent the collective emission from dust heated by nearby hot stars 
in star-forming regions, like, e.g. the clumpy component of the 


Popescu et al. (2000a 2011 1 models. Heating from star-forming 


regions in the galactic plane of NGC 891 was indeed found to be 
dominant in the shorter FIR bands ( Hughes et al.|2014l l. 

In Fig. ^ we show the scaleheights of the thin (black dia¬ 
monds) and thick (red squares) components obtained fitting ob¬ 
servations in region X with this three-component model and 
compared them with the scaleheights obtained using the stan¬ 
dard two-component model (black and red dots). The addition 
of the super-thin component leads to a milder gradient in the 
scaleheight of the thin component for wavelengths in the range 
70 - 250 pm. This has consequences on the dust temperature es¬ 
timation leading to a flatter temperature profile as a function of 


the vertical height above the midplane, more similarly to that 
expected from RT models. 

On the other hand, we find that the 8 and 24 pm data are 
best fit with a null contribution from the super-thin disc, thus not 
leading to any changes in the scaleheight of the thin and thick 
components. If we admit that 70 - 250pm data are affected by 
unresolved emission even in region X while 8 and 24 pm data are 
not (or to a low extent), then the estimate of the fraction of small 
grains destroyed at high galactic latitudes (~ 2.5 kpc above the 
galactic disc) given in Section 5.1 will represent a lower limit. 

However, the inclusion of a super-thin disc to the standard 
two-component vertical profile adds complexity to the model 
and leads to overfitting the data. In fact, the resolution of the 
available instruments is not sufficient to probe scales of tens of 
parsecs at the distance of this galaxy. We have therefore pre¬ 
sented the results from fitting a three-component model only 
qualitatively in order to show the consequences of the presence 
of unresolved sources in the midplane. 
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Fig. 8: Scaleheights of the thin (black diamonds) and thick (red 
squares) dust components as obtained from fitting observations 
in region X including the super-thin component. Green triangles 


represent the dust scaleheights as estimated by Xilouris et al. 


(1999i. Black and red dots are the scaleheights as retrieved in 
region X using our standard two-component model. 


6.2. Heating sources 


In the case where region X was not highly affected by unresolved 
emission the big difference in dust temperature between the thin 
and thick components requires two different heating sources. 
The thin component is likely heated by the galactic radiation 
field, which would drop off much more quickly than expected 
by radiative transfer models. On the other hand, the thick dust 
component might be heated by the presence of a hot halo gas or 
from strong radiation from nearby hot evolved stars. 

Assuming that the hot halo gas is distributed as observed 
bylHodges-Kluck & Bregman ( 2013| l and following the method 
by jBocchio et al~|^2013| l we computed the dust emission in the 
case of collisional heating only. We obtained a dust temperature 
Fdust = 19,18.2 and 17.4 K at z = 1,2 and 3 kpc from the mid¬ 
plane, respectively, which is not enough to explain the observed 
high dust temperature. 
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Another possibility comes from the strong radiation field 
from evolved stars present in the halo. In fact, accordingly to 
the model by |Flores-Fajardo et al. ( 2011| l, hot low-mass evolved 
stars are present in the halo of NGC 891 and dominate the ion¬ 
isation of the halo gas and could possibly heat nearby dust to 
high temperatures. A similar source of heating was considered 
by |Burgdorf et al.| ( |2007| l in order to explain the observed emis¬ 
sion from halo dust. This model is also consistent with HST ob¬ 
servations from |Ibata et ah] ( |2009| l, which show a thick stellar 
component with a vertical scaleheight Zs = (1.44 + 0.03) kpc. 

Furthermore, observations of the Hi gas distribution from 
Oosterloo et al. (2007|l show the presence of an Hi halo ex¬ 


tending out to about 14 kpc from the midplane in the northeast, 
southeast, and southwest quadrants while in the northwest quad¬ 
rant a flare extending out to ~ 22 kpc is observed. A few years 
later Yim et al. ( |2011| l, using Very Large Array observations de¬ 
tected a thin and a thick Hi components. They modelled the thick 
distribution with a Gaussian function with a FWHM of ~ 2 kpc 
which roughly translates into an exponential decay with a scale- 
height of ~ 1.4kpc. Moreover, this galaxy has been mapped in 
the 1=2-1 and 1-0 lines of '^CO with the IRAM 30-m telescope 
by IGarcia-Burillo et ^ ( 1992|l and in the J= 1 -0 line with the 
Nobeyama 45-m telescope by |Sofue & Nakai ( |1993| l. These two 
teams detected a ‘halo’ of molecular gas > 7 kpc in radius and 
2-3 kpc thick. This component is very weak in the 1-0 line and 
seems to be relatively hot and not too optically thick. This means 
that dust. Hi gas, molecular gas and stars follow a similar distri¬ 
bution in the halo and that dust might be embedded in atomic / 
molecular gas and heated by this secondary distribution of stars. 

On the contrary, if we admit that, even in region X there 
might be some contamination due to the presence of a super- 
thin disc (as suggested by our three-component model analysis 
in Sect |6.1| l, then the gradient in scaleheight of the thin com¬ 
ponent could be heavily smoothed. This would lead to a flatter 
profile in temperature, possibly similar to that of the thick com¬ 
ponent. This scenario then does not require two distinct sources 
of heating and the galactic radiation held as modeled by various 
radiative transfer models would be sufficient to explain the high 
dust temperature in the halo of the galaxy. 


6.3. The origin of extra-pianar dust 

We have detected dust emission at z > 2 kpc from the galac¬ 
tic disc. However, the evidence for extra-planar dust is not lim¬ 
ited to NGC 891. Dusty clouds have been observed in extinction 
up to z ~ 2 in the thick discs of the majority of the edge-on 
galaxies (|Howk & Savage|1997[ 1999[ 2000 Rossa et al.|2004[ 


[Thompson et al.||2004[ |Howk||2()05| |26l2| l. The optical images 

of these galaxies often show filamentary structures perpendic¬ 
ular to the galactic disc, linking the thin and thick discs. Also, 
extended (up to ~ 6 kpc) mid-infrared aromatic feature emission 
is seen around a number of galaxies (see e.g. Engelbracht et al. 


|2006[|Irwin & Madden|2006[|McCormick et al.|2013|l. 

The presence of dust in galactic haloes can be the result of 
several important phenomena. We identified at least three mech¬ 
anisms, possibly able to explain the presence of dust at these 
distances above the disc. 1) Dust can be expelled due to galac¬ 
tic fountains or stellar feedback (local winds). Observational and 
theoretical evidence show that small grains are partly destroyed 


during coupling with shocks and galactic fountains (Borkowski 


et al. 

2006 

IWilliams et al.|2006t|Micelotta et al.|20101|Bocchio 

et al. 

2014 

1 , in agreement with what we found. 2) Slow global 


winds (e.g. cosmic ray driven winds) could also explain the pres¬ 
ence of dust in the galactic halo. These winds are slowly and 


continuously emitted by spiral galaxies pushing dust grains to 
high galactic latitudes. Also, the expansion of these winds lead 
to a rapid drop in the gas temperature, therefore reducing grain 
sputtering ( jPopescu et al.|2000b I. 3) Another possibility comes 
from the radiation pressure. In fact, Ferrara et al. ( 1991| l show 
that this mechanism might displace dust up to several kpc from 
the galactic disc. 4) Finally, a significant contribution to the ob¬ 
served dust might come from dust formed in-situ around AGB 
stars present in the halo ( |Howk|2012| ). However, this scenario is 
not very likely since it would not explain the grain erosion above 
the galactic midplane. 

A significant amount of dust is observed even further out in 
the IGM around galaxies ( Zaritsky||1994[ Menard et al.||2010| l. 
The thick dusty disc described here would be an interface from 
which dust is ejected into the IGM. 


7. Summary 

We have analysed deep PACS data, together with IRAC, MIPS 
and SPIRE data and we have detected dust emission from the 
halo of NGC 891. Our main findings are summarised here. 


Dust emission at all wavelengths is best fit by a double ex¬ 
ponential profile revealing the presence of two dust compo¬ 
nents: a thin and a thick disc. 

The measured scaleheight of the thin disc shows a gradient 
through the FIR. This can be possibly explained by a dra¬ 
matic drop in dust temperature within the thin disc scale- 
height. Alternatively, contamination from a super-thin disc 
may also produce this effect and reconcile the observations 
with standard RT models for this galaxy. 

We obtain a precise measurement of the thick halo scale- 
height, Zd ,2 = (1.44 + 0.12) kpc, which is in agreement with 
previous estimates (jBurgdorf et al.||2007| [Kamphuis et al.j 


2007)[Whaley et al.|2009||Seon et al.|2014|l. 

We And that a non-negligible fraction (2 - 3.3%) of the mass 
of dust is located further than 2 kpc from the midplane, in 
accord with findings by [Seon et al. (2014 1 . Furthermore, 
the relative abundance of small grains with respect to large 
grains drops from imsg/imlo = (6.2 + 0.1) x 10“^ in the 
galactic disc to iwsg/wAg = (3-9 + 0.1) x 10“^ at a height 
of ~ 2.5 kpc above the plane suggesting that dust is hit by in¬ 
terstellar shocks or galactic fountains and entrained together 
with gas to a height of a few kpc from the midplane. 

The detected halo dust component is found to have a similar 
extent of the atomic / molecular gas distribution and of the 
thick stellar disc described by Ibata et al. ( 2009| l. The thick 
stellar disc can therefore be a possible source of heating. 
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Appendix A: A simple geometrical model 


For our analysis we considered strips of -2.5',-t-2.5' from the 
centre of the galaxy and parallel to its major axis. In this region 
the radial profile is relatively constant (emission always larger 
than ~ 1/10 of the peak) and we therefore considered the me¬ 
dian emission to compute the vertical profile of the galaxy. In 
this appendix we use a simple geometrical model to show the 
general validity of our approach, and study the bias in the esti¬ 
mates induced by geometrical components we did not account 
for in the fitting. 

As an example we consider PACS 100 fim data throughout 
this section. We construct an array with a pixel size equal to 
the corresponding instrument pixel size (i.e. 1") with flat null 
background. We perform a fit to the radial emission profile and 
assumed it as the radial profile of the intrinsic emission distri¬ 
bution. We use the intrinsic vertical profile extracted in region X 
(see Sect. 3.4 Z(j,i = 0.114kpc and Zd ,2 - 1.02kpc) as an esti¬ 


mate of the real vertical profile and we added it to all the pixels 
around the midplane. The simple model constructed in this way 
is then smoothed to the instrument resolution convolving it to 
the corresponding PSF (the B(mod) PSF for PACS 100 /tm) and 
adjusted in order to match the observed surface brightness. 


Table A.l: Scaleheights of the thin and thick components as 
extracted (at galactic scale and in region X) from observations 
(obs.) and from our geometrical model. Case a and case y6 refer 
to the model without and with the inclusion of a super-thin disc, 
respectively. 



galactic 

region 

X 

Zi,l 

(kpc) 

Zd,2 

(kpc) 

Zd.l 

(kpc) 

Zd,2 

(kpc) 

obs. 

0.087 ± 0.003 

0.63 ± 0.02 

0.114 ±0.004 

1.02 ±0.10 

a 

0.113 ±0.004 

0.98 ± 0.06 

0.111 ±0.004 

1.01 ±0.10 

p 

0.069 ± 0.002 

0.87 ± 0.06 

0.111 ±0.004 

1.04 ±0.07 


We extract vertical profiles by averaging over the whole 
galactic disc and in region X. In Fig. [A.l | we show the vertical 
profile extracted in region X from our model (cyan solid line) 
and compare it to the observed vertical profile from region X 
(black solid line). We then consider a two-component intrinsic 
exponential profile (Eq. [^, convolve it with the instrument LSF 
and fit it to the vertical profile of the convolved model (like¬ 
wise in Sect. |3.4[ ). The resulting profile is illustrated by the green 
dashed line in Fig. A. 1 and the parameters retrieved from the fit 
are reported as case a in Table [AT| The scaleheight of the thin 
and thick components for both regions are compatible with those 
used to construct our simple geometrical model. The compatibil¬ 
ity between these results demonstrates that the assumption of a 
flat radial profile does not introduce any systematic error on the 
estimation of the vertical scaleheights. 

However, this simple model is not able to reproduce the dis¬ 
crepancy in scaleheight fitting observations at galactic scale or 
considering only region X, which could be due to unresolved 
sources not accounted for in our fitting procedure (see Sect. 
and |6.1[ ). To test this hypothesis, we made two simple tests. 

In the first test, we include a few (5) isolated point sources 
along the galactic plane, to simulate the main features in the ob¬ 
served surface brightness distribution. These features could rep¬ 
resent inhomogeneities in the disc due to e.g. spiral arms seen 
in projection (see the discussion in [Bianchi & Xilouns]|2011||. 
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z (kpc) 

Fig. A.l: Normalised vertical profiles for PACS 100 /rm. Blue 
and black solid lines refer to the observed profiles at galactic 
scale and for region X, respectively. Cyan and orange solid lines 
represent the vertical profiles obtained from our simple geomet¬ 
rical model without (case a) and with (case j3) a super-thin disc 
included, respectively, while green and red dashed lines indicate 
the results of the fit in the two cases. 


Being point-like, their vertical profile would follow that of a PSF, 
and this would make the average profile decline faster than in 
the case of two diffuse discs. The intensity of the point sources 
(modelled as bright pixels) and of the diffuse discs where scaled 
to match the observed average profile, and an average profile 
was extracted from the model as we have done before. Though 
the new profile (not shown) is narrower than the profile for re¬ 
gion X, it is very close to it and far from the observed averaged 
profile. Thus, isolated point sources are not likely to be the main 
cause of the discrepancy. 

In a second test, we assume that half of the emission in the 
midplane comes from a very thin diffuse disc which is not re¬ 
solved vertically (and thus has a vertical profile similar to that 
of the LSF). This is analogous to dividing the intrinsic profile 
by a factor two for all the pixels which are not in the midplane. 
An exception is made for the pixels around region X where no 
bright source is added. We then convolve this model with the 
instrument PSF and adjust the normalisation to the observed 
surface brightness. The vertical profile at galactic scale of the 


convolved model is shown by the orange solid line in Fig. A. 1 


and compared to the observed vertical profile at galactic scale 
(blue solid line). With the same method used above, the scale- 
height of the thin and thick components are retrieved from fitting 
our convolved model. The resulting vertical profile is shown in 
Fig. |A.l|in red dashed line and fitting parameters are reported in 
Table A. 1 (case j3). The inclusion of an unresolved disc reduces 
the scaleheight of both the thin and thick components leading 
to a good agreement between the convolved and observed verti¬ 
cal profiles at galactic scale. Thus, the presence of a super-thin 
disc can explain the discrepancy between the vertical profiles 
extracted from observations at galactic scale and in region X. 

It has to be noted that the scaleheight of the thin and thick 
components retrieved from fitting our convolved model in region 
X are not affected by the inclusion of unresolved sources along 
the midplane (see Table |A.1[ case j3, region X). However, this 
does not imply that we can exclude that part of the emission 
from the midplane of region X is actually due to the presence of 
unresolved sources. In fact, this represents a scenario that we are 
not able to test with this simple geometrical model. 


Appendix B: Beyond “standard” Herschel PSFs. 

In this work we attempt to take into account the major effects 
that could modify the PACS and SPIRE PSFs. In particular, we 
study the variations induced by colour corrections (B.l), we de¬ 
rived a new set of “observed” PSFs for PACS (B.2) and finally 
tested the impact of these modifications on the derived vertical 
scaleheights (B.3). 

B.1. PSF modification according to the source spectrum 

The effective frequency, Veff, of an instrument equipped with 
a narrowband filter corresponds to the frequency at the centre 
of the filter. For larger filters, Veff depends on the spectrum of 
the observed source and therefore the instrument PSF shrinks 
or stretch according to the source spectrum. In the latter case, 
a PSF estimated from observations of an astronomical object is 
only valid if the source that we observe has a similar spectrum to 
that of the object used for the PSF characterisation. However, the 
spectrum of an asteroid or a planet can be very different from the 
dust emission spectrum. For this reason we modified the PACS 
and SPIRE PSEs according to the typical dust temperature and 
values. 

Eollowing the SPIRE Handbook (2014) a PSF depends on 
the frequency of the observed source as: 

B(0,v) = B(0x(—f ,Veff), (B.l) 

where 9 is the radial distance from the centre, v is the source 
frequency, Vefj is the effective frequency for which the PSF has 
been estimated and y is a parameter which takes into account the 
slope of the spectrum within a given band (y = -1 for an ideal 
telescope). 

The SPIRE ICC provides values of y and Veff for all SPIRE 
bands, in particular y = -0.85 and Veff = 1215 GHz for SPIRE 
250 pm. 



Eig. B.l: PACS 70, 100 and 160 pm (blue, green and red lines, 
respectively) and SPIRE 250 pm (black lines) PSEs, standard 
and modified assuming Tjust = 20K and yfijust = 2. 


Eor PACS bands we assume y = -1 and we estimate Vgff 
from the spectrum of Vesta and Mars as follows: 


vS(v)T(v)W(v)dv 
S(y)T(y)W(y)dy ’ 
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where S (v) is the source of the object used for the estimation of 
the PSF, T{v) is the transmission function and W(v) the bolome¬ 
ter respons^ 

Finally we compute the PSF for a given spectrum, S'(v) of 
the object that we want to observe (i.e. dust typically at F = 20 K 
and ~ 2): 


or non-linear behaviour of the central pixels, therefore allowing 
a good estimate of the core of the PSF. On the other hand, Mars 
is much brighter and leads to the complete saturation of the cen¬ 
tral pixels but gives much more information on the PSF faint 
wings. The combination of these observations can thus lead to 
an accurate estimate of the PSF. 


B(0, T,I3) 


T(v)S'(v)W(v)B(0,v)dv 
f T(y)S'(y)W(y)dy 


(B.3) 


For a more complete description of the method used the reader 
is referred to |Bocchio|(|20 1 41l . 

In Fig. |B.l|we show PACS 70, 100 and 160 jum and SPIRE 
250 jum PSFs for the standard case and modified according to 
the typical dust emission spectrum. Here, the “standard” PSFs 
for PACS are B (std) while for SPIRE are A (std). We note that 
for all wavelengths variations are never more relevant than 1% 
level. 

Given the lack of information on the IRAC and MIPS 
monochromatic PSFs and given the low impact of the colour 
corrections for PACS and SPIRE PSEs we do not perform this 
work for Spitzer PSEs. 


B.2. PACS PSFs estimation from observations 



Fig. B.3: PSF Radial profiles as estimated from observations 
(solid lines) and modeled by the ICC (dotted lines) for PACS 
70 pm (blue lines), 100 pm (green lines) and 160 pm (red lines). 



Fig. B.2: Radial profiles of Vesta (black solid line) and Mars 
(black dotted line) observations at 70 pm. In the inset the ratio 
between the Vesta and Mars profiles is shown (red solid line). 
Regions are indicated (see text for details). 


In order to obtain a good estimate of an instrument PSF it 
is necessary to have information both on its core and on faint 
structures in the outskirts. However, this involves a wide range 
in intensity and observations of a single object are often not suf¬ 
ficient to properly characterise a PSF. 

The PACS ICC provides an optical model for the instrument 
PSF taking into account aberrations of the real telescope (“as 
built” PSlQ. 

We use observations of an asteroid, Vesta, and a planet, Mars, 
to compare the observed PSF with that modeled by the ICC. 
Vesta is a sufficiently faint object not to lead to any saturation 

’ PACS and SPIRE filter transmission functions and bolometer re¬ 
sponse can be retrieved can be accessed from the HIRE environment. 

* http://herschel.esac.esa.int/twiki/pub/Public/ 
PacsCalibrationWeb/PACSPSF_PICC-ME-TN-029_v2.O.pd£ 


Observations of Vesta (ObsID: 1342195472, 1342195473 
for 70 pm, 1342195476, 1342195477 for 100 pm and 
1342195472, 1342195473, 1342195476, 1342195477 for 160 
pm) and Mars (ObsID: 1342231157 - 1342231160 for 70 
pm, 1342231161 - 1342231164 for 100 pm and 1342231157 
- 1342231164 for 160 pm) were reduced following the same 
method as described in Section|2] 

First, all the images have been rotated in a way to obtain 
the Z-axis of the telescope during observations pointing upward. 
Then, for all the images we removed an average background flux 
computed in regions sufficiently far from the source. We nor¬ 
malised the images of Vesta to unity and rescaled the images of 
Mars so to have the same flux in an intermediate region between 
the core and the outskirts of the PSF. 

The average radial profiles of Vesta and Mars for PACS 70 
pm are shown in Fig. |B.2| The inset illustrates the ratio between 
the flux of Vesta and Mars as a function of the distance from 
the PSF centre. We note in a region close to the centre (region 
1), the radial profile extracted from Mars does not follow that of 
Vesta because of saturation / non-linear issues, while at distances 
> 70" (region 3) the Vesta profile is not detected at a significant 
level. At intermediate distances (region2, from xi ~ 20" to X 2 ~ 
50") the ratio is close to unity. As an estimate of the PSF we use 
Vesta observations in region 1, Mars observations in region 3 and 
the average value between the two for region 2. We performed 
the same operation for PACS 100 and 160 pm, using xi - 30 
and X 2 - 50 for PACS 100 pm and xi - 40 and X 2 - 70 for 
PACS 160 pm. 

In order to obtain a smooth PSF in xi and X 2 we adopt the 
following smooth step function (Ebert et al.|2002|l: 


fix) = 6 


X - Xi 
X2 - Xi 


- 15 


X- Xi 
X2 - Xi 


+ 10 


X - Xl 
X2 - Xl 


(B.4) 


which has zero first and second order derivatives at the bound¬ 
aries Xl and X 2 . 
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In Fig. |B.3| we show radial profiles of the PACS PSFs as 
released by the ICC and obtained in this way. Differences are 
evident at all wavelengths and distances. In particular, ICC PSFs 
are systematically narrower by 10-30% than those we obtained 
from observations. 

B.3. Scaleheights assuming different PSFs 

With the method described in Section [3 we estimated the scale- 
height of the dust emission profile at different wavelengths. The 
instrument PSF plays a key role in this calculation. In Fig. |B.4| 
we show the scaleheights as retrieved from the fitting proce¬ 
dure using the different PSFs we considered in the main text. 
We note that the effect of modifying the PSF according to the 
observed source spectrum has no relevant effect on the estima¬ 
tion of the scaleheight, reflecting the small effect on the PSF (see 
Section [BTT] i. On the contrary, passing from the “as built” PACS 
PSF (A) to that estimated from observations (B) importantly af¬ 
fects the dust scaleheight. 
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Fig. B.4: Scaleheights of the thin (black symbols) and thick 
(red symbols) dust components assuming different PSFs. Green 


et al. 


squares represent the dust scaleheights as estimated by Xilouris 

799 ^ . 
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